param c{LINKS} default 1;

var P {FLOWS, LINKS} binary default 0;
var w {LINKS} default 0;

# min c^t w
minimize w_cost:  (sum {(i1,i2) in LINKS} c[i1,i2]*w[i1,i2]) /2 ;

# P B^T = D
s.t. mass_baW {(r1,r2) in FLOWS, n1 in NODES}:
  sum {(i1,i2) in LINKS} P[r1,r2,i1,i2]*B[n1,i1,i2] 
      = D[r1,r2,n1];

s.t. cap_aggrW {(i1,i2) in LINKS}:
  w[i1,i2] = sum{(r1,r2) in FLOWS} P[r1,r2,i1,i2]*M[r1,r2];

s.t. work_sym {(r1,r2) in FLOWS, (i1,i2) in LINKS}:
  P[r1,r2,i1,i2] = P[r2,r1,i2,i1];
# P is symatric for symetric flows

problem find_work:  w_cost,  mass_baW, cap_aggrW, work_sym, P, w;
